Hi Luca,

Thanks for this PR! I'll let Steven judge the astrophysics, but code quality looks fine to me overall. I do have two more software-technical points.

First, you added tests, which is great! However, they plot their results, which means that a human would have to look at those plots to see if they make sense. We'd really like to be able to run the tests automatically (there are a bunch of tests in src/amuse/tests actually that are set up for this), for which we need some assert statements in the test that check whether the result is as expected, and which crash otherwise.

Would it be possible to add those?

I'm currently nearing the end of a big project to replace the AMUSE build system, and that also reorganises those tests in src/amuse/tests, so it's fine to have them as separate scripts in Tests/. I'll merge them in with the rest when that new build system lands so that they can be run automatically with everything else, ensuring that your additions won't be broken by some future change.

Second, the fix to the loop causes some repetition, and I think this part could be written more compactly, but I need to have a more careful look at a more reasonable hour. I'll post if I find a better solution.

Hi Lourens,

Thank you for your comments!

Concerning the update to the evolve function, I have not found a way to make it more compact than the proposed version while ensuring that the next timestep is not set to an unintended small value. But if you find a more compact way of doing it, that’d be great!

Concerning the tests, I understand your point. I've had a look at the test_mesa_15140.py script in src/amuse/test/suite/codes_tests which contains some tests with assert statements. I have added two similar tests at the end of this script (test26 and test27).

Test 26 checks the updated evolve_until method, in particular that after several calls to evolve_for, the next timestep has the expected value (i.e. 1.2* the value given as argument to evolve_for).

Test 27 checks the functionality of the added getters. For this one I have used the already written Test4, and instead of getting the luminosity, mass and temperature (as done in Test4), I use the added getters.

Should I do an other PR to include the updated version of test_mesa_15140.py ?

Hi Luca,

If you make a new commit with the new tests on your main_update_mesa_r15140_interface branch, and push it to GitHub, then they will show up here automatically and be included in this PR. So no need to make a new one.

Hi Lourens,

I have just made a new commit which deletes the two tests I had added in the original PR, and modifies the test script test_mesa_15140.py script in src/amuse/test/suite/codes_tests. The new version of the script contains the two additional tests.

Great! I've had a look at that loop in the meantime, and I think it can be rewritten along these lines:

# pseudocode
while star_age < end_time - epsilon:
    dt_save = dt_next
    if star_age + dt_next / secyer < end_time - epsilon:
        dt_next = (end_time - star_age) * secyer
    do_evolve_one_step()
    dt_next = dt_save

That avoids the need for the separate last step, reducing duplication. Do you think that will work?

Thank you for this pseudo-code.

I however think it does not provide the intended result. Let me try to explain why with an example. In some case, the call to evolve_until will be done with a pretty large value for the argument delta_t (i.e. substantially bigger than the current MESA timestep, so that MESA would need several timesteps to reach end_time. The name evolve_until is actually quite ambiguous, since we don't evolve until delta_t but for delta_t.

Example: star_age = 1 Gyr, call to evolve_until with argument delta_t = 1 Gyr.
end_time = star_age + delta_t = 2 Gyr.

But the timestep required by MESA could be much lower, for instance ~1 Myr, so that several steps are required to reach end_time.

1. With the code you proposed and in my example, we enter the loop since star_age < end_time
2. dt is saved
3. star_age + dt_next = 1 Gyr + 1 Myr < 2 Gyr so we enter the if
4. But here, dt_next is directly set to end_time - star_age i.e. 1 Gyr
5. We evolve one step and then directly reach end_time (since dt_next has been defined this way).

So in any case instead of a loop there is only one iteration before reaching the final age. We can not simply set MESA timestep this way, the code can note accomodate for such important modification (in particular increase) of the timestep. So I don't think this works.

In the way I wrote the loop, it lets MESA choose independently the value it wants for the timestep as long as the chosen value does not make the star reach an age higher than end_time. Only for the last step it sets the timestep to reach exactly end_time. It then sets dt_next back to the value it was the step before, so that if evolve_until is called again, the first timestep is what it was just before the last step of the previous loop.

There will be situations were the loop is not even entered, when dt_next (timestep wanted by mesa )> delta_t (timestep given as argument to evolve_until). In this case, the loop is not entered, and the timestep given to mesa is just delta_t.

Tell me if these explanations help. There may however be a way to make the code more compact keeping the same results, I can try to find a way to make it more compact if you want.

Ah, I got the test in the if the wrong way around! I meant to set dt_next only if we're going to step past the end. So what I intended to write was

while star_age < end_time - epsilon:
    dt_save = dt_next
    if end_time - epsilon < star_age + dt_next / secyer:
        dt_next = (end_time - star_age) * secyer
    do_evolve_one_step()
    dt_next = dt_save

I didn't realise that do_evolve_one_step() could update dt_next, but it's obvious in retrospect, why have it be part of the state otherwise? Would this work instead?

while star_age < end_time - epsilon:
    if end_time - epsilon < star_age + dt_next / secyer:
        dt_save = dt_next
        dt_saved = true
        dt_next = (end_time - star_age) * secyer
    do_evolve_one_step()
    if dt_saved:
        dt_next = dt_save

with an additional dt_saved logical variable? Alternatively, dt_save could start out as -1 and we could check if it's not negative to see if it's been set. Or actually, the dt_next = dt_save bit could go after the loop and always run, assuming that that epsilon guarantees that we always run the saving operation on the last loop iteration?

Thank you for these suggestions!

Following these lines I have updated the evolve_until (last commit) accordingly. I have essentially written it as you suggested. The only difference with what you suggested is that actually dt_next is systematically set at the last step so that star_age reaches exactly end_time (and therefore dt_save is also set). So the condition to set dt_next to dt_save at the end of a loop should not be that dt_save has been set since it happens in any case at the end of a loop, but that there was more than one call to evolve_one_step.

Essentialy the if statement will always be true at the end of a loop, so dt_save will systematically be set. The condition to put back dt_next to dt_save at the end of a loop should be that at least one time the if statement was false (i.e. there was more than one call to evolve_one_step).

This is because if only one call to evolve_one_step was required to reach end_time, the MESA timestep has not been reduced more than necessary to reach end_time (it could not have a larger value since at least one step was required to reach end_time) and in this case after the loop dt_next already has the "correct value" it should have, i.e. 1.2*previous_time_step.

I hope these explanations make sense.

I have thus modified the code to keep this idea, but it is now written in a more compact way following your suggestion. I have tested this version with the test I provided in the previous commit and it works fine.